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FLUID: A NUMERICAL INTERPOLATION PROCEDURE FOR OBTAINING 
THERMODYNAMIC AND TRANSPORT PROPERTIES OF FLUIDS 
by Theodore E. Fessler 
Lewis Research Center 

SUMMARY 

A computer -program subroutine, FLUID, has been developed to calculate the ther- 
modynamic and transport properties of pure fluid sub stances. It provides for deter- 
mining the thermodynamic state from assigned values for temperature -density, 
pressure -density, temperature-pressure, pressure -entropy, or pressure -enthalpy. 
FLUID considers liquid or two -phase (liquid-gas) conditions as well as the gas phase. A 
van der Waals model is used to obtain approximate state values; these values are then 
corrected for real -gas effects by model -correction factors obtained from tables based on 
experimental data. Saturation conditions, specific heat, entropy, and enthalpy data are 
included in the tables for each gas. Since these tables are external to the FLUID sub- 
routine itself, FLUID can implement any gas for which a set of tables has been gen- 
erated. (A setup phase is used to establish pointers dynamically to the tables for a spe- 
cific gas. ) Data-table preparation is described. FLUID is available in both SFTRAN 
and FORTRAN. 


INTRODUCTION 

For several years the computer code GASP (ref. 1) has been used at the Lewis Re- 
search Center to supply thermodynamic and transport property values for engineering 
calculations. Because current and projected use of GASP is so great, the replacement 
code FLUID has been developed as a cost-saving measure. FLUID provides thermody- 
namic and transport properties in good agreement with GASP and WASP and some sav- 
ings in computer time. 

GASP is a FORTRAN subprogram, the current version of which has been imple- 
mented for 10 gases: hydrogen, oxygen, nitrogen, argon, carbon dioxide, carbon mon- 
oxide, fluorine, helium, neon, and methane. A similar program, WASP (ref. 2), is 
used for water or steam properties. In GASP, a modified virial equation of state is used 


to express pressure as a function of temperature and density. In WASP, thermodynamic 
properties are calculated from the Helmholtz free-energy function. In both programs , 
transport properties are calculated from empirical functions in temperature and density. 
GASP and WASP are applicable to both the liquid and gas phases, including saturation 
conditions in which a mixture of phases occurs. 

All of the thermodynamic properties can be obtained from an equation of state (or 
from the Helmholtz function in the case of WASP). Hence, in principle, all of the ther- 
modynamic property data for one gas can be combined by a constrained, weighted least- 
squares technique in which the equation -of -state coefficients are obtained numerically. 
The coefficients used in GASP (and WASP) were obtained in this way. 

A disadvantage to this approach is that it is hard to find a functional form for the 
equation of state (or the Helmholtz function) that works well over the entire temperature - 
density plane. Many terms are required. And, it is costly to solve these functions 
inversely - for independent variables in terms of dependent variables. For example, 
GASP takes much longer to compute temperature and density from pressure and enthalpy 
than it does to compute pressure and enthalpy from temperature and density. Another 
problem may occur if the choice of the algorithm used to obtain these inverse solutions 
depends on the argument values. If that be the case, not only will the program be more 
complex, but there may be discontinuities in calculated values for the different 
algorithms. 

In the FLUID subroutine, using a simple gas model and empirical corrections to it, 
obtained by numerical interpolation, avoids some of these complexities and makes pos- 
sible some savings in computer time. 


SYMBOLS 


A sonic velocity 

a tabulated argument values 

C specific heat 

C j , . . . , C 8 constants used in data tables 
c coefficients used in interpolation formulas 

d reduced density, p/p 

F,G tabular functions 

H enthalpy 

P pressure 

p reduced pressure, P/P 

L 
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R 
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T 
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U 
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W 

x 
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gas constant 

entropy 

temperature 

reduced temperature, T/T c 
internal energy 
specific volume 

mass fraction of one phase in a two-phase system 
argument of interpolation 
model correction factors 
constants 

ratio of specific heats, Cp/C y 
density 


Subscripts: 

c 

g 

l 

P 

t 

v 

0 


value at critical point 

gas -phase value 

liquid-phase value 

at constant pressure 

at constant temperature 

at constant volume 

function of temperature for ideal gas 


PROGRAM DESIGN 

FLUID, like GASP, is designed to calculate thermodynamic and transport properties 
of pure fluids in both the liquid and gas phases. In FLUID, an equation of state provides 
a model for calculating thermodynamic properties; tables of model-correction factors 
are used to bring the equation of state into agreement with accepted values, as for exam- 
ple, those from GASP. Other tables similar to the model -correction -factor tables are 
used directly, without a model, to calculate viscosity and thermal conductivity. Inter- 
polation formulas are used for argument values between those in the tables. 
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Two auxiliary subroutines are used by the FLUID subroutine: NTRP, in which all 
interpolations are performed, and CUBIC, in which the equation of state is solved for 
density. 

The tables for each gas are external to FLUID. New gases may be implemented by 
building the appropriate tables from an already existing code like GASP or by some di- 
rect use of tabulated values. The tables of thermodynamic properties are in nondimen- 
sional form to allow for a change in the system of units. 


Equation of State 

The van der Waals equation of state can qualitatively represent the behavior of a 
real fluid in both its liquid and gaseous states. Hence, it has obvious application in the 
work at hand. 

Model . - In its reduced form, van der Waals’ equation is 

(p + 3d 2 )(3 - d) = j 
8td 

(See SYMBOLS. ) Real liquids and gases show a behavior somewhat different from that 
predicted by van der Waals' equation. We choose to account for these differences by us- 
ing a correction factor Zj, where 

(p . + ifj 2 H 3 ~ d) = Z. (t,d) 

8td 1 

which is found experimentally to be a smooth, slowly varying function of temperature and 
density. (Similar methods have been used by others to obtain correction factors, resid- 
ual values, or excess functions that relate observed property values to values predicted 
by a model. See, e. g. , ref. 3. ) If values of Z^ are available on a coarse tabular grid 
in t and d, two-dimensional interpolation can be used to obtain Z^ for intermediate 
values of t and d. 

The behavior of real fluids will also depart from van der Waals' model whenever two- 
phase conditions occur. In that case, pressure, liquid -phase density, and gas-phase 
density are expressible as functions of temperature only. These functions can be ar- 
ranged in tabular form so that interpolation can supply function values at intermediate 
argument values. 

Limitations of model. - Although the van der Waals equation gives a satisfactory 
qualitative representation of the behavior of a real gas, it fails in one important aspect 
for the present application. This failure comes from the derivation of the model, in 
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which it is assumed that molecules have some hard, nonzero volume. That is, there is 
a certain density that cannot be exceeded: three times the critical density. But some 
real fluids have liquid densities that are greater than this model maximum. For exam- 
ple, the density of water at room temperature is 3.15 times the critical density. This 
limits the liquid -phase range of FLUID to temperatures greater than 400 K for water (al- 
though gas and two -phase calculations can be performed down to 273 K). 

One way around this limitation is to assume a somewhat higher, pseudocritical den- 
sity, both when the model correction tables are generated and when they are used. For 
water, a pseudocritical density of 0. 34 g/cm instead of the observed value of 0. 317 
could be used. This has been found to work well, except that the resulting tables are not 
to be used near the critical point. In this example for water, this would be the region in 
which 350 K < T < 410 K and 0. 25 g/cm 2 < p < 0. 40 g/cm 2 . Fortunately, many fluids 
of interest are not so seriously affected by this density limitation as is water. 


Other Thermodynamic Functions 


The van der Waals equation of state implies certain behavior of other thermodynamic 
variables and fluid properties as well. Since the specific heat at constant volume 
satisfies 


1 

T 




it follows that, for a van der Waals gas, C y must be at most a function of temperature: 


C v “ C v,0^ 


Also, since 


C 


P 




it can be shown that 



R 

d(3 - d) 2 
4t 
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Then, since entropy can be expressed as 


it follows that 



dV + Constant 


S = S Q (t) + R In 



where SpCt) is a function of temperature only for an ideal gas. Likewise, internal en- 
ergy can be expressed as 


U 




- P 


dV + Constant 


which leads to 


U = U 0 (t)-iRT c d 


And, from this result and H = U + PV, we can obtain 


H = U n (t) + — RT f — - d] 


8 


,3d 


Finally, the sonic velocity in a van der Waals' gas can be derived from 



to be 


A = 



8t 

3-d 



3-d 



where y = Cp/C v . 

These equations only approximate the behavior of real liquids and gases. As in the 
case of the equation of state, we choose to account for these differences by the use of 
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correction factors, which can be put in tabular form. The following relationships 
define these correction factors: 


S = Z 2 (t,d)x 


S Q (t) + R In (L 


H = Z 3 (t,d) x 


u ° (t) + 5 RT ‘(i H 


C v = Z 4 (t,d)xC v>0 (t) 


c 


p 


Z 5 (t,d)x 


C v + 


R 

1 d 0 - d) 2 
4t - 


A = Z 6 (t,d) 



If saturation occurs, two-phase conditions prevail and we can express liquid- and 
gas -phase values of entropy and enthalpy as functions of temperature only. (We do not 
define specific heats or sonic velocity when two phases occur. ) The entropy and enthalpy 
of the mixture are then given by 


S = W 7 S 7 + W S 
11 g g 


H = W 7 H 7 + W H 
Z l g g 

where the subscripts l and g refer to the liquid and gas phases, respectively, and W 
is a mass fraction in one phase. 


Method of Interpolation 

Interpolation is used in this package to obtain values of functions and correction fac- 
tors at points intermediate to tabulated values. The method used has been called the 
"osculatory" or ' 'double -Lagrangian ’ ' interpolation method. Following Lagrange's three- 
point interpolation formula, we find the value of a function by 
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F(x) = 


(x - a 0 )(x - a,) (x - a,)(x - aJ (x - a.)(x - a 0 ) 

— F( ai ) + — F(a 2 ) + * L_ F(a„) 

( a ^ - a 2 ^ a l - a 3 ^ (&2 ~ a i ^ a 2 "" a 3 ^ ( a 3 - a ]^( a 3 ~ a 2 ) 

where x is an argument intermediate between the tabular values a^, a 2 , and a 3 and 
F(a^), etc., are the tabulated values of the function on the a grid. Suppose x lies be- 
tween a 2 and a 3 . Then F(x) could also be calculated by using tabular values at a^, 
a 3 , and a^(if a^ exists). Let the formula using a^, a^, and a 3 be called LOWER(x) 
and the formula using a 2 , a 3 , and a^ be called UPPER(x). Then by double interpola- 
tion, 

G(x) = — — — LOWER(x) + — — UPPER(x) 

(32 - a 3 ) (a 3 - ag) 

in which a weighted average of the two interpolations is prescribed. This method of in- 
terpolation results in a very smooth function because G(x) has a continuous first 
derivative. 

The interpolation formulas just described can be generalized as 

n 

G(x) = 2 c- x F(a.) 
i=l 1 1 

where n = 3 for the three-point formula and n = 4 for the weighted-average formula. 
The values of c^ depend on where x falls in the table of arguments but not on the func- 
tion values in the table. Thus, if several tabular functions share the same argument 
table, the same interpolation coefficients apply equally to all functions for any given 
value of x. In order to take advantage of this, calls to the interpolation subprogram 
NTRP are made in two steps: first, the coefficients c- are calculated and saved and 
then, for each tabular function to be interpolated, a call is made with that function table 
as the argument. If the argument value x falls in the argument table so that a 2 < x 
< a * , where m is the length of the argument table, the weighted-average formula can 
be used. If not, the three-point formula is used. 

For interpolation in two dimensions, the interpolation formula can be expressed as 

o 

a double summation over n tabulated values. In this case, each argument table has 
associated with it a set of c- values that are calculated by using these formulas. 
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Solutions for Other Independent Variables 


Three sets of tabular data are used in this package: (1) the ideal -gas functions 
s 0 (t), U Q (t), and C v , Q (t); (2) saturation values p(t), dj(t), d g (t), S z (t), S g (t), H z (t), and 
H g (t); and (3) correction factors Z^t, d). Interpolation in these tables provides all data 
necessary to calculate thermodynamic properties for a given temperature and density. 

But more generally, since any two thermodynamic variables determine the fluid state, 
we require solutions for those cases where the given variables may be other than t and 
d. Four of these cases have been implemented and will now be considered individually. 

Pressure and density. - When density is an input parameter, a preliminary calcula- 
tion must be made to determine if two-phase conditions prevail. If the pressure is less 
than critical, the liquid- and gas -phase densities can be obtained from the saturation 
tables by interpolation. A two -phase condition is indicated when d g < d < d^ . If that is 
the case, the liquid- and gas -phase values of the other thermodynamic variables can also 
be obtained by interpolation of the saturation tables. 

If only one phase exists, the temperature is found by successive approximation, with 
Z^ = 1. 0 as a starting value. Iteration continues until temperature differences are 
smaller than a convergence criterion, Cg. 

Temperature and pressure . - If the temperature is less than critical, van der Waals' 
equation yields three roots when solved for density. The smallest of these corresponds 
to the gas value, and the largest corresponds to the liquid value. (These values are cal- 
culated in subroutine CUBIC. ) If the pressure is less than saturation, the gas value is 
used; otherwise, the liquid value is used. Van der Waals’ equation yields only one root 
(gas) when the temperature is greater than critical. 

In either case, density is found by successive approximation, with Z^ = 1. 0 as a 
starting value. Iteration continues until density differences are smaller than a conver- 
gence criterion, Cg. 

Pressure and entropy. - When entropy is an input parameter, a preliminary calcula- 
tion is made to decide if two-phase conditions prevail. If the pressure is less than crit- 
ical, the liquid and gas entropies can be obtained from the saturation tables by interpola- 
tion. A two -phase condition is indicated when < S < S g . If that is the case, the 
liquid- and gas -phase values of the other thermodynamic variables can also be obtained 
by interpolation of the saturation tables. 

If only one phase exists, density and temperature are found by successive approx- 
imation. If the given pressure is less than critical, the liquid or gas saturation density 
is used as a starting value, depending on whether the entropy is less than or more than 
critical. If the pressure is more than critical, the calculation starts at critical density. 
By using the procedure described for pressure and density, trial values of density lead 
to trial values of temperature. 
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The entropy calculated from trial values of temperature and density is then com- 
pared with the given entropy, and a new trial value of density is computed by using an 
approximation for entropy as a function of density. If the fluid state is gas -like (entropy 
greater than critical), this approximation is of the form S(d) = a + /3 log(d), where a 
and /3 are constants; if the fluid state is liquid -like, the approximation is of the form 
S(d) = a + /3 log (3 - d). Iteration continues until entropy differences satisfy a conver- 
gence criterion, C^. 

Pressure and enthalpy . - The method used here is similar to that for the case of 
pressure and entropy; a preliminary calculation detects two -phase conditions and, if 
that is the case, the other thermodynamic variables are obtained by interpolation of the 
saturation tables. 

If only one phase exists, successive approximation is used to solve for temperature 
and density. The fluid state is considered gas-like if the given enthalpy is more than the 
saturated -liquid enthalpy (and the pressure is less than critical) or if the given enthalpy 
is more than the critical -point enthalpy (and the pressure is greater than critical). The 
approximations for enthalpy as a function of density are of the form H(d) = a + /3/d if 
gas-like and H(d) = a + /3d if liquid-like. Iteration continues until enthalpy differences 
satisfy a convergence criterion, Cg. 


Calling Conventions for FLUID 

FLUID, NTRP, and CUBIC were coded in SFTRAN (Structured FO RTRAN ), a pro- 
gramming language recently introduced at the Lewis Research Center (ref. 4). This 
language was chosen for its self -documenting features. It allows and even encourages 
the grouping of instructions into small, natural units that can be given unique, descrip- 
tive names. In listings, these are displayed in a manner that makes their internal struc- 
ture immediately apparent to the eye. 

Other than the manner in which branching and looping are handled, SFTRAN is very 

much like FORTRAN; SFTRAN routines can call, or be called by, FORTRAN routines. 

The SFTRAN language has been implemented as a precompiler that generates FORTRAN 

source code from SFTRAN source code. The FORTRAN versions of FLUID, NTRP, and 

CUBIC so generated have been included in the FLUID package for use by those not having 

an SFTRAN precompiler. (This FORTRAN code will not be helpful in understanding 

FLUID; the FORTRAN produced by the SFTRAN precompiler is hard to read. ) 

# 

Initializing FLUID for a new gas . - FLUID does not contain any data tables. Instead, 
all of the tables for one gas are collected into a gas -tables subroutine (one subroutine for 
each gas). Then, when the user's program is ready to calculate properties for a par- 
ticular gas, it calls the gas -tables subprogram for that gas, which in turn calls FLUIDO, 
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setting up pointers to the appropriate tables. (See appendix A for a description of the 
entries in FLUID. ) Subsequent calls to FLUID for property values will then use these 
tables until such time as a different gas-tables subprogram is called by the user’s 
program. 

A library of gas-tables subroutines has been developed. Each subroutine in the li- 
brary has been named for the particular gas it implements. Hence, 02 is for oxygen, 
PH2 is for para -hydrogen, etc. These subroutines are available in source code and can 
be transported to other computers. 

Calls to FLUID for property values. - Five combinations of input variables have 
been implemented: (T,p), (P,p), (T,P), (P, S), and (P,H). A particular one is specified 
by integer ENTRY. Temperature, pressure, and density are transmitted explicitly as 
calling arguments. Entropy and enthalpy are transmitted in the array PROPS, along with 
the specific heats, sonic velocity, and the transport properties. An integer, NP, is used 
to specify the length of PROPS and hence how many of these properties are to be 
calculated. 

As an example, consider that at some point in the user 's program, a calculation of 
oxygen properties is required. In particular, suppose it is desired to obtain the temper- 
ature, density, and specific heats at constant volume from pressure and enthalpy. For 
this, the calling program might contain the sequence 
CALL 02 
NPROPS=4 
ENTRY=5 
PROPS(3)=H 

CALL FLUID (T, P, D, PROPS, NPROPS, ENTRY, VAPOR, ERROR) 

IF (ERROR. NE. 0) DO (REPORT ERROR IN FLUID) 

IF (VAPOR) THEN 

DO -(SPECIAL CASE: MIXED-PHASE CONDITION) 

ELSE 

S=PR0PS(2) 

CV=PR0PS(4) 

DO (NORMAL-CASE CALCULATION) 

END 

(It is assumed in this example that ENTRY and ERROR have been specified as INTEGER 
type and VAPOR as LOGICAL type and that PROPS is a dimensioned array. ) In this ex- 
ample, the highest index in PROPS of a desired property is 4 (specific heat at constant 
volume), and so the number of properties requested from FLUID need be only 4. The 
remaining properties (5 to 8) would not be calculated, thereby effecting a saving of com- 
puter time. The reader is referred to appendix B for more examples of the use of 
FLUID. 

If the values of the input variables correspond to a mixed-phase condition, VAPOR 
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will be. TRUE. In that case, specific heats, sonic velocity, and transport properties 
are not calculated but liquid- and gas -phase values of mass fraction, density, entropy, 
and enthalpy may be available in the COMMON block/FLUIDC/depending on how many 
properties (NP) have been asked for. Liquid- or gas-phase values of specific heats, 
sonic velocity, and transport properties can be obtained by the auxiliary entries FLUIDL 
and FLUIDG. 

Error reports . - FLUID recognizes certain kinds of errors in calling arguments. 
Fatal errors (due to an error in coding in the calling program) are announced on I/O 
UNIT 6 and result in a program STOP. Errors due to one or more of the table arguments 
being out of range are not fatal; they are signaled by certain bits in integer ERROR. If 
for some reason, convergence is not achieved in solving for temperature and/or density 
from other input variables, this condition is also signaled by a bit in ERROR. 


GENERATION OF TABLES 

If all, or nearly all, of the property data for a particular fluid are available in a 
program such as GASP, making tables for FLUID can be a programmed operation. On 
the other hand, it is possible (but with a greater effort) to make the desired tables di- 
rectly from tabular data. In the example used here, tables for oxygen were produced 
from data partly from GASP and partly from NBS Circular 564 (ref. 5). 


Gas Constant, Critical Values, and Convergence Constants 
The following values were used for the oxygen tables: 


R (C x ) = 0.25983 J/g-K 
T c (C 2 ) = 154. 78 K 
P c (C 3 ) = 5. 083 MN/m 2 
p c (C 4 ) = 0.4325 g/cm 3 


c 5 = 0. 001 

C g = 0. 002 

c 7 = 0. 001 

Cg = 0. 001 


The dimensions used for R, T, P, and p also determine the dimensions used for entropy, 
enthalpy, specific heats, and sonic velocity. For the dimensions indicated, entropy and 
specific heats will be in J/g-K, enthalpy will be in J/g, and sonic velocity will be in J/g 
to the 1/2 power. (If the more usual units of m/sec are desired for sonic velocity, the 
user should multiply the value returned from FLUID by 158. 33. ) The four convergence 
constants required by FLUID (Cg to Cg) were determined by trial. Their values are not 
critical. The particular values chosen give good accuracy and at the same time min- 
imize the number of iterations required to obtain that accuracy. 
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Saturation Tables 


FLUID does not make use of a theoretical vapor -pressure equation such as the 
Clausius-Clapeyron equation. Instead, it uses the interpolation formulas described ear- 
lier. Since these formulas are cubic in the argument, they do not approximate well the 
vapor -pressure and vapor-density data. But if temperature arguments for the saturation 
tables are chosen to make each vapor -pressure value about twice that of the preceding 
value, the interpolated vapor pressures will track the correct function within about 1 part 
per thousand. (The interpolation algorithm requires temperature arguments in all tables 
to be in ascending order, but they need not have equal intervals. This requirement is 
imposed by the table -lookup procedure used in subroutine NTRP. ) The reduced- 
temperature argument values used for oxygen are 0.42, 0.43, 0.45, 0.475, 0.5, 0.53, 
0.56 , 0.60 , 0.65 , 0.70 , 0.75 , 0 . 80 , 0 . 85 , 0.90 , 0.95 , 0.97 , 0.982 , 0.99 , 0.9 95, 

0.998, 0.999, and 1.0. 

Inasmuch as entropy and enthalpy are defined as integrals, a reference value may be 
assigned arbitrarily to each. But for FLUID, it is necessary that entropy and enthalpy 
be positive in the temperature and density range covered by the model-correction tables. 
Otherwise, the use of model -correction factors to obtain real-gas values of entropy and 
enthalpy will not work. This requirement is satisfied by using reference values that 
make the saturated-vapor entropy and enthalpy twice their saturated -liquid values at the 
lowest temperature in the saturation tables. 

The last value of each saturation-table quantity is the critical value. The critical 
values of reduced temperature, density, and pressure are unity, of course; the critical 
entropy and critical enthalpy are calculated as the arithmetic average of the next -to -last 
liquid- and gas -phase values. 


Tables in Temperature Only 


FLUID requires tables of values for the functions SQ(t), UQ(t), and C y Q(t), together 
with their argument values. According to the procedure discussed in the section Other 
Thermodynamic Functions, these ideal -gas functions are evaluated for some low -density 
value by 


S Q (t) =S(t,d) -Rln 
U 0 (t) = H(t,d)-|RT c ^-d) 
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Cv.ofO'CyM) 


in which the S, H, and C y values are obtained from the data base. (The density chosen 
should be low enough that saturation conditions are avoided over the entire temperature 
range used. ) 

The values of the argument t used for these tables should be selected to give good 
agreement between interpolated values and measured values. If specific heat shows 
marked variations with temperature, it may be necessary to use finer spacing where that 
occurs. The reduced -temperature arguments used for oxygen are 0.42 , 0. 5 , 0. 6, 0. 72, 
0.85, 0.98, 1. 1, 1.25, 1.5, 1.75, 2.0, 2.5, 3.0, 3.7, 4.75, 6.0, 7.2, and 8.0. 

The version of GASP available at the time of this writing implements oxygen data up 
to 500 K (reduced temperature of 3. 23). It was desired to extend FLUID'S capability for 
oxygen to a reduced temperature of 8. 0 (1238 K). To this end, the C y Q(t) table was ex- 
tended by using the 0. 01-atmosphere specific heats for oxygen of reference 5 above 
300 K. Then the SQ(t) values were calculated by numerical integration of 


S 0 (t 2 ,d) 

R 


SoM* 

+ 

R 




dt 


Rt 


and the Up(t) values were calculated by numerical integration of 

Po(t2 . 1 ) . ¥l' a ) , f 2 C v, 0 ^ , / p c 3\ p<'2> - ay 

RT c RT c 1 R U T c R 8 /L d J 

r l 


(For an ideal van der Waals' gas, the last term on the right is zero; for oxygen, it was 
cal ciliated by using the van der Waals model and the low -density Z j value for oxygen of 
1.2832.) 


Tables in Temperature and Density 

The t-d tables make up the bulk of the data used by FLUID (92 percent in the case 
of oxygen). Some care is needed in selecting the values for temperature and density if 
table size is to be minimized while maintaining good agreement with the data base over 
the entire temperature -density plane. The values used for oxygen are given in table I. 
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TABLE I. - REDUCED -TEMPERATURE 


AND REDUCED -DENSITY ARGUMENT 
VALUES USED IN TWO-ARGUMENT 


TABLES FOR OXYGEN 


Reduced temperature, 
t 

Reduced density, 
d 

0. 42 

0.0001 

.45 

.0003 

. 5 

.0007 

.59 

.002 

.68 

.01 

. 8 

.03 

.92 

. 1 

.985 

.3 

1. 05 

.6 

1. 15 

1.0 

1.3 

1.45 

1.6 

1. 85 

2.0 

2.2 

2.6 

2.46 

3.5 

2.65 

4.6 

2.79 

6.0 

2.87 

8.0 

2.92 


These values were arrived at by a trial-and -error process in which some general 
guides were developed. For instance, it was found best to have adjacent argument- 
intervals in density differ by at most a factor of 3. Also, it was found helpful to have 
some temperature -density points distributed more or less uniformly alongside the sat- 
uration curve. Thus, at t = 0.42, the points at d = 0. 0003 (gas) and d = 2. 92 (liquid) 
lie just outside the two-phase region; at t = 0. 5, the points at d = 0. 0007 and d = 2. 87 
lie just outside the two-phase region, etc. For other gases, circumstances may well 
require a different mesh; table I should be used only as a starter. 

Model -correction factors . - Even though GASP was not implemented for oxygen at 
temperatures above 500 K, it produced Z. values that were in agreement with refer- 
ence 5 up to 1238 K. (This is not too surprising; none of the Z. values varies more 
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than about 0. 002 from 500 K to 1238 K for pressures below 30 MN/m^. ) Hence, Z^ 
values for the entire table were based on values calculated by using GASP. At some high- 
density, high- temperature points, the GASP values were obviously erroneous. (These 
points were all beyond the pressure range of the data on which GASP is based. ) In these 
few cases, reasonable Zj values were filled in just to make the tables complete. 

The Z^ can be calculated by using the defining equations given previously. Values 
are obtained at each mesh point in the grid of temperature -density arguments except 
those for which saturation occurs. FLUID does not attempt to calculate the Z^ for sat- 
uration conditions, but the interpolation algorithm requires that values be there. For 
this purpose, the Z^ tables under the saturation line can be filled in for each density 
argument by extrapolating downward in temperature from the unsaturated region. 

Viscosity and conductivity . - In addition to the model -correction factors, the t-d 
tables can also be used for transport properties. In the oxygen tables, viscosity and con- 
ductivity values have been included. For these, no model is assumed. Instead, the two- 
dimensional interpolation in t-d gives directly the desired property. Values of vis- 
cosity and conductivity obtained from GASP were used for FLUID’S tables. The units of 
viscosity are g/cm-sec and the units of conductivity are W/cm-K. 


PROGRAM PERFORMANCE AND SIZE 

Central Processing Unit (CPU) times required to complete one calculation with 
FLUID depend on the entry type and on the number of additional properties called for. 
These times have been measured on the Lewis Research Center's IBM 360 computer sys- 
tem as (typically) 4 milliseconds for T,p entries, 16 milliseconds for P,p entries, 

12 milliseconds for T,P entries, and 50 milliseconds for P, S and P,H entries. 
Additional properties (NP > 1) require 1 millisecond each. These times are nearly in- 
dependent of data-table size. The CPU time required by an initialization call to FLUID 
is insignificant. 

One of the goals of the implementation of FLUID was to minimize program size. 

The Lewis Research Center’s IBM 360 computer is a virtual -memory, time-share sys- 
tem. In this kind of system, small program size makes more CPU time available in 
real (connect) time. The compiled code for the FLUID subroutine and the two subrou- 
tines it references, NTRP and CUBIC, require 4141 (4-byte) words of storage. 

The compiled code for the oxygen tables, developed in the manner described here 
and implemented as a subroutine, requires a total of 2960 words of storage. These 
tables provide the real-gas equation of state: the thermodynamic properties S, H, C y , 
Cp; and sonic velocity, viscosity, and thermal conductivity. The tables are sufficiently 
fine to provide an agreement with the data base, which is, for the most part, better than 
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1 part per thousand. Tables for other gases may be expected to fall within a factor of 2 
of this number of storage words for a comparable accuracy. 

Appendix B lists a sample program that demonstrates the use of FLUID with the 
oxygen tables. It can be used to obtain thermodynamic properties in units other than 
those for which the data tables for FLUID were originally prepared. For example, by 
using the constants built into the program, values will be obtained that may be compared 
directly with the oxygen values in reference 5. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, May 5, 1977, 

505-01. 
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APPENDIX A 


ABSTRACT OF ENTRY POINTS IN FLUID 


SUBROUTINE FLUIDO (G ,NG ,GG ,TT ,DD ,Nl ,N2 ,N3 ,SAT ,NS ,C ,L) 
INTITIALIZING ENTRY FOR FLUID-PROPERTIES PACKAGE. SETS UP TABLE 
POINTERS FOR A PARTICULAR FLUID. SUBSEQUENT CALLS (TO FLUID) USE 
THESE POINTERS. 


C CALLING ARGUMENTS: 

C G = (NON-DIMENSIONAL) TABLES IN TEMPERATURE ONLY. 

C G (1 , 1) = ARGUMENT VALUES (TEMPERATURE). 

C G (1 ,2) = ENTROPY, SO. 

C G ( 1 3 ) = ENERGY, UO . 

C G (1 , 4) = SPECIFIC HEAT AT CONSTANT VOLUME, CVO . 

C NG = LENGTH OF EACH G TABLE. 

C GG = TABLES IN TEMPERATURE AND DENSITY. 

C GG (1,1,1) = Zl, THE MODEL-DEPARTURE FACTOR. 

C GG (1,1,2) = Z2 FACTOR FOR ENTROPY (S) . 

C GG (1,1,3) = Z3 FACTOR FOR ENTHALPY (H) . 

C GG (1,1,4) = Z4 FACTOR FOR SPECIFIC HEAT (CV) . 

C GG (1,1,5) = Z5 FACTOR FOR SPECIFIC HEAT (CP). 

C GG (1,1,6) = Z6 FACTOR FOR VELOCITY OF SOUND. 

C GG ( 1 , 1 ,N) = OTHER FACTORS (EG. TRANSPORT PROPERTIES). 

C TT = TEMPERATURE ARGUMENTS (REDUCED) FOR GG TABLES. 

C DD = DENSITY ARGUMENTS (REDUCED) FOR GG TABLES. 

C Nl = LENGTH OF TT TABLE AND FIRST DIMENSION OF GG TABLE. 

C N2 = LENGTH OF DD TABLE AND SECOND DIMENSION OF GG TABLE. 

C N3 = THIRD DIMENSION (NUMBER OF) GG TABLES. 

C SAT = (NON-DIMENSIONAL) SATURATION TABLES. 

C SAT (1,1) = SATURATION TEMPERATURE. 

C SAT (1,2) = SATURATION PRESSURE. 

C SAT (1,3) = SATURATED LIQUID DENSITY. 

C SAT (1,4) = SATURATED GAS DENSITY. 

C SAT (1,5) = SATURATED LIQUID ENTROPY. 

C SAT (1,6) = SATURATED GAS ENTROPY. 

C SAT (1,7) = SATURATED LIQUID ENTHALPY. 

C SAT (1,8) = SATURATED GAS ENTHALPY. 

C NS = LENGTH OF EACH SATURATION TABLE. 

C C = MISCELLANEOUS CONSTANTS: 

C C (1) = SPECIFIC GAS CONSTANT, R. 

C C (2 ) = CRITICAL TEMPERATURE, TC . 

C C (3) = CRITICAL PRESSURE, PC. 

C C (4) = CRITICAL DENSITY, DC (OR PSEUDO VALUE FOR BEST RANGE) . 

C C (5) = CONVERGENCE REQUIREMENT FOR TEMPERATURE. 

C C (6) = CONVERGENCE REQUIREMENT FOR DENSITY. 

C C (7 ) = CONVERGENCE REQUIREMENT FOR ENTROPY. 

C C ( 8 ) = CONVERGENCE REQUIREMENT FOR ENTHALPY. 

C L = MEMORY (4 WORDS) FOR TABLE INDICES. 

C VARIABLES IN COMMON BLOCK /FLUIDC/: 

C GAMMA = RATIO OF SPECIFIC HEATS, CP/CV. 

C WL = MASS-FRACTION IN THE LIQUID PHASE. 

C WG = MASS-FRACTION IN THE GAS PHASE. 

C DENSL = DENSITY OF SATURATED LIQUID, SAME UNITS AS DC. 

C DENSG = DENSITY OF SATURATED GAS, SAME UNITS AS DC. 

C ENTL = ENTROPY OF SATURATED LIQUID, SAME UNITS AS R. 
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C ENTG = ENTROPY OF SATURATED GAS, SAME UNITS AS R. 

C ENTHL = ENTHALPY OF SATURATED LIQUID, SAME UNITS AS R*TC . 

C ENTHG = ENTHALPY OF SATURATED GAS, SAME UNITS AS R*TC . 


COMMON /FLUIDC/ GAMMA ,WL ,WG ,DENSL ,DENSG ,ENTL , ENTG , ENTHL , ENTHG 


ENTRY FLUID (TEMP , PRES , DENS , PROPS ,NP , ENTRY , VAPOR, ERROR) 

C MAIN ENTRY FOR PROPERTY CALCULATIONS. 

C CALLING ARGUMENTS: 

C TEMP = FLUID TEMPERATURE, SAME UNITS AS TC . 

C PRES = FLUID PRESSURE, SAME UNITS AS PC. 

C DENS = FLUID DENSITY, SAME UNITS AS DC. 

C PROPS = OTHER FLUID PROPERTIES: 

C PROPS (1) = COMPRESSIBILITY, PRES/ (R*DENS*TEMP) . 

C PROPS (2) = ENTROPY, SAME UNITS AS R. 

C PROPS (3) = ENTHALPY, SAME UNITS AS R*TC. 

C PROPS (4) = SPECIFIC HEAT, CV, SAME UNITS AS R. 

C PROPS (5) = SPECIFIC HEAT, CP, SAME UNITS AS R. 

C PROPS (6) = SONIC VELOCITY, SAME UNITS AS SQRT (R*TC) . 

C NP = NUMBER OF PROPERTIES TO BE CALCULATED. 

C ENTRY = INTEGER THAT SPECIFIES WHICH VARIABLES ARE INPUT: 

C = 1 IF TEMPERATURE AND DENSITY ARE GIVEN. 

C = 2 IF PRESSURE AND DENSITY ARE GIVEN. 

C = 3 IF TEMPERATURE AND PRESSURE ARE GIVEN. 

C = 4 IF PRESSURE AND ENTROPY ARE GIVEN. 

C = 5 IF PRESSURE AND ENTHALPY ARE GIVEN. 

C VAPOR = .TRUE. IF THE FLUID IS SATURATED. IN THAT CASE, 

C VALUES OF LIQUID AND GAS PHASES ARE GIVEN IN THE 

C COMMON BLOCK /FLUIDC/. 

C ERROR = ERROR FLAGS (BITS — LEAST SIGNIFICANT = 1) : 

C ** IF ERROR = 0, ALL IS WELL ** 

C BIT 1 = OUT OF RANGE IN SAT TABLE. 

C BIT 2 = OUT OF RANGE IN G TABLE. 

C BIT 3 = OUT OF RANGE IN TT TABLE. 

C BIT 4 = OUT OF RANGE IN DD TABLE. 

C BIT 5 = CONVERGENCE NOT ACHIEVED IN SOLVING INVERSE FUNCTION. 


ENTRY FLUIDL (PROPS ,NP , ERROR) 

C ENTRY FOR GETTING LIQUID-PHASE PROPERTIES NOT IN COMMON. 

ENTRY FLUIDG (PROPS ,NP , ERROR) 

C ENTRY FOR GETTING GAS -PHASE PROPERTIES NOT IN COMMON. 
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APPENDIX B 


SOURCE -CODE LISTING (IN SFTRAN) OF A TEST PROGRAM FOR 
OBTAINING OXYGEN VALUES FROM FLUID PACKAGE 


C TESTS FLUID SUBPROGRAM AND OXYGEN TABLES. 

INCLUDE (TYPE STATEMENTS, ETC.) 

INCLUDE (CONSTANTS FOR COMPARING OXYGEN RESULTS TO NBS 564) 

NAMELIST /INPUT/ TEMP , PRES , DENS , ENT ,ENTH , ENTRY , REF ,NP , TO ,P0 , 

* TCONV , PCONV , DCONV , SCONV ,HCONV ,CCONV , ACONV , 

* S REF , HREF 


C MAIN PROGRAM FLOW. 

CALL 02 
DO WITH 

READ (5 , INPUT, DONE=END) 

UNTIL (DONE) 

IF ( ENTRY. GE.l .AND. ENTRY. LE. 5) THEN 
DO CASE (ENTRY, 5) 

CASE 1 

DO (TEMP AND DENS INPUT) 

CASE 2 

DO (PRES AND DENS INPUT) 

CASE 3 

DO (TEMP AND PRES INPUT) 

CASE 4 

DO (PRES AND ENT INPUT) 

CASE 5 

DO (PRES AND ENTH INPUT) 

END 

IF (VAPOR) THEN 

WRITE (6,1) TEMP, PRES, DENS, ENT, ENTH 
ELSE 

WRITE (6,2) TEMP, PRES, DENS, ENT, ENTH ,CV, CP, A 
IF (NP.GT.6) WRITE (6,3) (PROPS (I) ,1=7 ,NP) 

END 

ELSE 

IF (VAPOR) THEN 

IF (ENTRY. EQ. 6) THEN 

DO (LIQUID-PHASE VALUES) 

ELSE 

IF (ENTRY. EQ. 7) THEN 

DO (GAS -PHASE VALUES) 

ELSE 

WRITE (6,4) 

END 

END 

ELSE 

WRITE (6,4) 

END 

END 

END 

STOP 
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DEFINITION (CONSTANTS FOR COMPARING OXYGEN RESULTS TO NBS 564) 


FOR TEMPERATURE IN K, PRESSURE IN ATMOSPHERES, REDUCED DENSITY, 
ENTROPY AS S/R, ENTHALPY AS H/(R*T0), SPECIFIC HEATS AS CV/R AND 
CP/R, AND REDUCED SONIC VELOCITY. (THE VALUES OF SREF AND HREF 
HERE GIVE AGREEMENT BETWEEN FLUID AND NBS 564 AT 273 K, 1 ATM.) 

DATA TO ,P0 ,TCONV ,PCONV / 0.0, 0.0, 1.0, 9.8692 /, 

DCONV , SCONV , HCONV , CCONV / 699.8, 3.8487, 0.01409, 3.8487 /, 
ACONV, SREF, HREF, NP / 0.10045, 1.180, 408.14, 6 / 

END 


C PROCEDURES . 

PROCEDURE (TEMP AND DENS INPUT) 
TSAVE=TEMP 
DSAVE=DENS 

DO (GET ALL PROPERTIES VIA FLUID) 

TEMP=TSAVE 

DENS=DSAVE 

END 

PROCEDURE (PRES AND DENS INPUT) 
PSAVE=PRES 
DSAVE=DENS 

DO (GET ALL PROPERTIES VIA FLUID) 

PRES=PSAVE 

DEN S=D SAVE 

END 

PROCEDURE (TEMP AND PRES INPUT) 
TSAVE=TEMP 
PSAVE=PRES 

DO (GET ALL PROPERTIES VIA FLUID) 

TEMP=TSAVE 

PRES=PSAVE 

END 

PROCEDURE (PRES AND ENT INPUT) 

PSAVE=PRES 

SSAVE=ENT 

DO (GET ALL PROPERTIES VIA FLUID) 

PRES=PSAVE 

ENT=SSAVE 

END 

PROCEDURE (PRES AND ENTH INPUT) 
PSAVE=PRES 
HSAVE=ENTH 

DO (GET ALL PROPERTIES VIA FLUID) 

PRES=PSAVE 

ENTH=HSAVE 

END 

PROCEDURE (GET ALL PROPERTIES VIA FLUID) 
T= (TEMP+T0) /TCONV 
P= (PRES+P0) /PCONV 
D=DENS/DCONV 
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H <n ro in \o r- 


S=ENT/ SCONV + SREF 
H=ENTH/HCONV + HREF 

CALL FLUID (T ,P ,D , PROPS ,NP , ENTRY , VAPOR , ERROR) 
IF (ERROR. NE.O) WRITE (6,5) ERROR 
IF (REF) THEN 
SREF=S 
HREF=H 
REF=. FALSE. 

END 

TEMP=T*TCONV - TO 
PRES=P*PCONV - PO 
DENS=D*DCONV 
ENT= (S-SREF) *SCONV 
ENTH= (H-HREF) *HCONV 
IF (VAPOR) THEN 
DENSL=DL*DCONV 
DENSG=DG*DCONV 
ENTL= (SL-SREF) *SCONV 
ENTG= (SG-SREF ) *SCONV 
ENTHL= (HL-HREF) *HCONV 
ENTHG= (HG-HREF) *HCONV 
ELSE 

C V=C V * CCON V 
CP=CP*CCONV 
A=A*ACONV 

END 

END 

PROCEDURE (LIQUID-PHASE VALUES) 

CALL FLUIDL (PROPS ,NP , ERROR) 

DO (REPORT PHASE VALUES) 

END 

PROCEDURE (GAS -PHASE VALUES) 

CALL FLUIDG (PROPS ,NP , ERROR) 

DO (REPORT PHASE VALUES) 

END 

PROCEDURE (REPORT PHASE VALUES) 

IF (ERROR. NE.O) WRITE (6,5) ERROR 

CV=CV*CCONV 

CP=CP*CCONV 

A=A*ACONV 

IF (ENTRY. EQ. 6) THEN 

WRITE (6,6) DENSL ,ENTL ,ENTHL ,CV ,CP ,A 
ELSE 

WRITE (6,7) DENSG , ENTG , ENTHG , CV , CP , A 

END 

IF (NP.GT.6) WRITE (6,3) (PROPS (I ) , 1=7 ,NP) 

END 

FORMAT (' T,P,D,S,H: ' , 5 ( 1PE13 . 4 ) ) 

FORMAT (' T,P,D,S,H,CV,CP,A: ' ,8 (1PE13 .4) ) 

FORMAT (' ADDITIONAL PROPS : ' , 8 (1PE13 . 4 ) ) 

FORMAT ( ' *ERROR* BAD VALUE FOR ENTRY ' ) 

FORMAT (' ERROR =', 13) 

FORMAT (' DL ,SL ,HL ,CV ,CP , A : ' , 6 (1PE13 .4 ) ) 

FORMAT (' DG,SG,HG,CV,CP,A: ' , 6 (1PE13 .4) ) 
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DEFINITION (TYPE STATEMENTS, ETC.) 

INTEGER ENTRY, ERROR 
LOGICAL VAPOR, REF, DONE 
DIMENSION PROPS (6) 

COMMON /FLUIDC/ GAMMA ,WL ,WG ,DL ,DG , SL , SG ,HL ,HG 
EQUIVALENCE (Z , PROPS(l)), (S , PROPS(2)), (H, PROPS(3)), 
(CV, PROPS (4) ) , (CP, PROPS (5) ) , (A, PROPS (6)) 

END 

END 
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